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ABSTRACT 

We show that diffusion due to chaotic mixing in the Neighbourhood of the Sun may 
not be as relevant as previously suggested in erasing phase space signatures of past 
Galactic accretion events. For this purpose, we analyse Solar Neighbourhood like vol¬ 
umes extracted from cosmological simulations that naturally account for chaotic or¬ 
bital behaviour induced by the strongly triaxial and cuspy shape of the resulting dark 
matter haloes, among other factors. In the approximation of an analytical static triax¬ 
ial model, our results show that a large fraction of stellar halo particles in such local 
volumes have chaos onset times (i.e., the timescale at which stars commonly associ¬ 
ated with chaotic orbits will exhibit their chaotic behaviour) significantly larger than 
a Hubble time. Furthermore, particles that do present a chaotic behaviour within a 
Hubble time do not exhibit significant diffusion in phase space. 

Key words: chaos: galaxies - galaxies: dynamics - methods: variational chaos indi¬ 
cators - methods: N- body simulations 


1 INTRODUCTION 

In galactic dynamics, the term chaos refers to the expo¬ 
nential divergence of initially nearby orbits in phase space. 
In non-integrable systems, initially nearby stars in strong 
chaotic regions drift away from each other very quickly, 
thus losing memory of their common origin (e.g. 


Merritt 


fc Vallun||1996| |Contopoulos|[2002] |Efthymiopoulos et al. 


2008 


and references therein). Understanding whether this 
physical process plays a major role in shaping the stellar 
phase space distribution of the Solar Neighbourhood is of 
key importance, as it is within this distribution that we 
hope to gain a significant insight into the formation his¬ 
tory of the Galaxy (Helmi fe de Zeeuw|2000 Johnston et al. 


2008 Gomez et al.[20101. 

The characterisation of the formation history of our 
own Galaxy is a very ambitious undertaking by modern as¬ 


tronomy (Freeman & Bland-Hawthorn 20021. Understand¬ 


ing how the Milky Way evolved to become the galaxy we 
currently inhabit would allow us not only to explore our 


E-mail: nmaffione@fcaglp. unlp.edu. ar 


origins, but also to understand galaxy formation in a more 
general context (Helmi 2008). More precisely, it can allow us 
to test the current paradigm of galaxy formation and evolu¬ 
tion. This theory predicts that the present-day population 
of galaxies grew in mass by merging with smaller compan¬ 
ions. As their potential wells grew deeper, galaxies contin¬ 
ued to accrete gas that cooled, collapsed into a disk, and 


gave rise to most of their stellar component (White & Rees 


1978). In addition to this in-situ population, every galaxy is 


predicted to have a minor fraction of its stellar content as¬ 
sociated with merger events (Searle & Zinn 19781. The tidal 
force that a satellite experiences as it orbits its host may be 
strong enough to disrupt it significantly. As a result of this 
interaction, initially spatially coherent and extended stellar 
streams are formed (e.g. Ibata et al. |[l994| |200l) |Bullock| 


fc Johnston|2005| |Belokurov et al.|2006[|2007 ]Cooper et al. 


2010). These streams are thus fossil signatures of this forma¬ 


tion process and their identification is key to reconstructing 
the merger history of our Galaxy. 

Stellar streams associated with the most ancient accre¬ 
tion events are expected to populate the inner Galactic re¬ 


gions, in particular the Solar Neighbourhood (Helmi et al. 
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2003 Gomez et al. 20131. Unfortunately, dynamical times 


in these regions are relatively short and streams tend to mix 
rapidly, losing their spatial coherence (Helmi fe White|1999 


Helmi||2008 I. Even though galaxies are essentially collision 


less (Binney & Tremaine 1987), and thus streams become 


more clustered in velocity as they diffuse in configuration 
space, statistically significant phase space overdensities are 
needed to identify such streams in the Solar Neighbourhood. 
Clearly, the longer it takes for a stream to diffuse in configu¬ 
ration space, the larger the chances of identifying it in phase 
space. The rate at which a stellar stream spatially dissolves 
depends not only on the orbit of its progenitor satellite, but 


also on the properties of the host galactic potential (Helmi 
fe White||1999 Helmi |[2008 Vogelsberger et al.||2008 i. In a 

galactic potential where a stream is on a regular orbit, its 
local density decreases in time as a power law, thus rela¬ 
tively slowly. Yet the current paradigm of galaxy formation 
predicts strongly triaxial and cuspy dark matter haloes (see 
e.g. |Jing fe Suto||2002[ | Allgood et al.|[2006] |Vera-Ciro et al.| 


2011). These two characteristics are known to be a signifi¬ 


cant source of chaos in a galactic potential (Schawrzschild 
1993} Merritt & Fridman 1996| Siopis & Kandrup 2000[ 


Voglis et al. ||2002| |Kandrup fc Siopis||2003| |Kalapotharakos| 

et al.|2004 Muzzio et al.|2005 l. On a chaotic orbit, the local 

stream density decreases at an exponential rate. Further¬ 
more, chaotic orbits can significantly drift in the space of 
quantities that are otherwise approximately conserved, such 
as angular momentum (e.g. Poveda et al.|[l992 Schuster & 


Allen 1997 Valluri et al. 20131, thus hindering its identifi¬ 


cation at the present day. 

One of the main goals of ongoing and past astrometric, 
photometric and spectroscopic missions is to map the stel¬ 
lar phase space distribution in the Solar Neighbourhood. In 
addition to Gaia (Perryman et al. 2001 Lindegren et al. 


2008), which will measure positions and velocities of more 


than a thousand million stars, several other projects have 
provided and will continue to provide complementary infer¬ 


mation (e.g. LAMOST, Zhao et al.[2006 

Zwitter et al.|2008 

HERMES, Wylie-de Boer & Freeman!2010; APOGEE, Ma- 

jewski et al. 2010 

Gilmore et al. 2012 DESI, Levi et al. 

2013 

Takada et al. 

2014 

). A meaningful interpretation of 


tile degree of substructure found in the Solar Neighbour¬ 
hood requires a deep understanding of the role that chaotic 
mixing plays in shaping the underlying substructure’s phase 
space distribution. 

Since the identification and characterisation of chaotic 
orbits is a fundamental step towards this goal, efficient and 
accurate tools for this purpose are essential. A seminal con¬ 
tribution to the field of chaos detection was made by Lya¬ 
punov ( Lyapunov| 18921 when he introduced the idea behind 
the so-called Lyapunov Characteristic Exponents (LCEs). 
LCEs are theoretical quantities that provide a measure of 
the rate of local exponential divergence of two initially 
nearby orbits in phase space. Thus, the LCEs are a very 
convenient way to distinguish between chaotic and regular 
motion and, particularly, to characterise chaos. Of particular 
importance is the largest LCE (1LCE), which is the LCE in 
the direction for which these two orbits diverge most rapidly. 
Theoretically, the characterisation of an orbit according to 
its 1LCE is done based on its asymptotic behaviour at infin¬ 
ity. The Lyapunov Indicator (LI), on the other hand, refers 
to the finite-time version of the 1LCE. A numerical value of 


the LI very close to zero indicates regular behaviour whereas 
any positive value indicates chaotic motion. The inverse of 
the LI provides a measure of the timescale for the manifes¬ 
tation of the exponential divergence. In practice, numerical 
finite-time techniques based on the concept of local expo¬ 
nential divergence like the LI (see e.g. Benettin et al.|[l980 


Skokos 2010) are commonly considered to be chaos indica¬ 


tors. 

Nowadays, there are many chaos indicators in the liter¬ 
ature that were developed based on the idea of the LCEs. 
Among the most used and tested indicators we find the Fast 


Lega] 

to 

o 

o 

o 

Lega & Froeschle 2001; Guzzo et al. 2002; Lega 

et al. 

o 

t-H 

o 

CM 

. The reliability shown by the FLI in previous 


works (e.g. Maffione et al.|2011a|b Darriba et al.|20l2 Maf- 

fione et al.|2013 1 makes this chaos indicator an ideal tool to 

characterise the role of chaos in Solar Neighbourhood-like 
volumes. 

In this work we will take advantage of a variant of 
the FLI, the so-called Orthogonal Fast Lyapunov Indicator 
(OFLI; Fouchard et alJ20021 to evaluate the importance of 
chaos and chaotic mixing in shaping the phase space distri¬ 
bution of halo stars in the Solar Neighbourhood. 

The paper is organised as follows: we describe the mod¬ 
els and the techniques in Section [2] and include a short 
but comprehensive characterisation of the behaviours of the 
OFLI in Section [3] In Section [4] we present our results on 
the actual impact of chaos in erasing local signatures of ac¬ 
cretion events on the stellar halo phase space and revise the 
widespread assumption that chaos must inevitably lead to 
diffusion. We discuss and summarise our results in Section [5} 
Finally, in the Appendix we include a brief introduction to 
the concept of local exponential divergence, the definition of 
the LI and a more formal description of the phenomenon of 
stickiness and sticky orbits. 


2 METHODOLOGY 

In this Section we introduce the simulations and numer¬ 
ical tools used to characterise the chaotic nature of the 
stellar halo phase space distribution enclosed within Solar 
Neighbourhood-like volumes. 


2.1 Simulations 

We use five high-resolution, fully cosmological IV-body sim¬ 
ulations of the formation of Milky Way (MW)-like dark mat¬ 
ter (DM) haloes carried out using the parallel Tree-PM code 


GADGET-3 (an upgraded version of GADGET-2, Springel 
|2005[ ) by the Aquarius Project ( |Springel et al.||2008a|b| ~ 
Each halo was first identified withi n a large cosmological 
periodic box of 100 h -1 Mpc a side (Boylan-Kolchin et al. 
[2009} and then re-simulated using a multi mass particle 
zoom-in technique. These DM-only simulations were per¬ 
formed using the following cosmological parameters: matter 
(dark and baryon) density, = 0.25; dark energy density, 
S7 a = 0.75; normalisation of the power matter spectrum, 
ag = 0.9; scalar spectral index, n s = 1 and Hubble con¬ 
stant, Hq = 100 h km s _1 Mpc^ 1 = 73 km s -1 Mpc -1 , con¬ 


sistent with WMAP 1 - and 5-year constraints (Spergel et al. 
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2003 Komatsu et al.|[2009 l. The DM haloes for the Aquar¬ 
ius Project were selected to have masses ~ 1O 12 M0, com¬ 
parable to the MW, and to be relatively isolated at z = 0. 
They were identified using a Friends-of-Friends (Davis et 
al.|1985 1 algorithm and self-bound subhaloes were identified 


with SUBFIND (Springel et al. 2001). Each MW-like DM 


halo was re-simulated at a series of progressively higher res¬ 
olutions. The experiments presented in this work are based 
on the simulations with the second highest resolution avail¬ 
able, and their main properties are presented in Table [l] For 
a more detailed description of the simulations, we refer the 
reader to Springel et al. (2008a b). 

To model the formation and present day properties of 
their galactic stellar haloes, these simulations were post- 
processed with a semi-analytic model of galaxy formation 


(Cooper et al. 2010). This semi-analytic model consists of 


a set of coupled differential equations describing the evolu¬ 
tion of baryons and derives the mass accretion histories and 
phase space information of halo stars from the underlying 
A-body DM-only simulations. Processes such as star forma¬ 
tion, AGN feedback, stellar winds and chemical enrichment 
are introduced in the model through differential equations 
that are controlled via a set of adjustable input parameters. 
The parameters were set to simultaneously match a range 
of observable quantities such as the galaxy luminosity func¬ 
tions in multiple wavebands ( Baugh et al.|20 05 Bower et al. 
2006 Font et ak|2008 1. 


The approach followed by Cooper et al. (20101 uses the 


technique known as particle tagging. The idea behind this 
technique is to assume that the most strongly bound DM 
particles in progenitor satellites can be used to trace the 
phase space distribution of their stars. At every time-step 
the 1% most-bound DM particles were selected to trace any 
newly formed stellar population in each galaxy in the sim¬ 
ulation. This fraction was set such that properties of the 
satellite population at 2 = 0 are consistent with those ob¬ 
served for the MW and M31 satellites. As a result of this 
procedure, each tagged particle has a different final stel¬ 
lar mass associated to it. From now on we will refer to the 
tagged particles as stellar particles. The main properties of 
the resulting stellar haloes are summarised in Table [l] Note 
that our stellar halo masses (A/*) also include the mass as¬ 
signed to the bulge component in |Cooper et al. ( |2010| ); see 
|Gomez et al.| ( f2013| hereinafter G13) and references therein 
for further details. 


Cooper et al. (20101 showed that particle tagging meth¬ 


ods applied to simulations with sufficient resolution can be 
used to generate MW-like stellar haloes that reproduce var¬ 
ious observables regarding the structure and characteris¬ 
tics of the Galactic stellar halo and its satellite population. 
Furthermore, as shown by G13, for most simulated stellar 
haloes, the measured velocity ellipsoids at 8 kpc are in good 
agree ment with the estima te for the local Galactic stellar 
halcp| (JChib aTfc Beers| 2000). Nonetheless, it should be kept 
in mind that the dynamical evolution of the baryonic com¬ 
ponents of galaxies in these simulations is much simplffied. 
As previously discussed in G13, this likely has an effect on, 


1 The velocity ellipsoids at 8 kpc from the galactic centre were 
measured along the direction of the major axis of the DM—halo 
to increase the particle resolution. 


e.g. the efficiency of satellite mass-loss due to tidal stripping, 
or the satellite’s internal structural changes due to adiabatic 
contraction, and possibly even on its final radial distribution 
(Libeskind et al.|2010 Romano-Diaz 2011 Schewtschenko & 


Maccio 2011 Geen et al.|2013| ). Recently, Bailin et al. 1 2014 ) 
compared the stellar haloes formed in fully Smoothed Par¬ 
ticle Hydrodynamics (SPH) simulations of galaxy formation 
with DM-only simulations of the same initial conditions. 
They found that the resulting stellar haloes have different 
concentrations and internal structure due to the different 
kinematics that DM particles show with respect to their 
SPH counterparts. However, in the particle tagging scheme 


used by Bailin et al. (2014) only one tagging operation is 


performed per satellite, at the time of its infall to the main 


halo, whereas Cooper et al. (20101 tag stars continuously, as 
they are formed. The Cooper et al. (2010) approach permits 


stars to diffuse in phase space within their parent satellite 
before its disruption, and thereby reproduces SPH results 
more closely than the tagging scheme tested by |Bailin et al.| 
120141 (Le Bret et al. submitted). 


2.2 The galactic potential 

The computation of chaos indicators (hereinafter CIs) to 
study the dynamics of our MW-like stellar haloes requires 
the integration of the equations of motion coupled with 
the first variational equations. The latter are used to track 
in time the evolution of the separation between initially 
nearby orbits in phase space (see Appendix [Aj. Due to 
the very high resolution of our A-body simulations, using 
frozen representations of the underlying galactic potential 
(see for instance Valluri et ak|2013) becomes computation¬ 
ally expensive. Methods to approximate the underlying po¬ 
tential based on series expansions can be extremely accurate 


(e-g. 

Clutton-Brock| 1973 

Hernquist & Ostriker 1992; 

Wein- 

berg 

1999 Kalapotharakos et al. 2008; Lowing et al. 2011 

Vasiliev 

2013 

Meiron et al. 

2014 

i. However, a very large 


number of expansion terms are needed in order to justify 


this approach. For instance, in Lowing et al. (2011), the au¬ 


thors used a halo expansion method to accurately fit the 
potential of the halo Aq-A2. Their analysis showed that a 
force accuracy of less than 1% can be achieved using a se¬ 
ries expansion that includes all moments up to n,l = 20. 
The resulting potential contains 8000 terms, rendering the 
derivation of the first variational equations unfeasible. In¬ 
stead, we chose to approximate each galactic potential with 
a suitable analytic model. Note that in this work we are 
dealing with pure DM simulations. Thus, the dynamics of 
the stellar particles are governed by the potential of the DM 
halo. 

It is important to mention that, in this work, the A- 
body simulations are mainly used to extract the parame¬ 
ters of the underlying triaxial potentials and to sample the 
phase space distribution of Solar Neighbourhood-like vol¬ 
umes (see Section [2.3[ ). In other words, our purpose is not 
to accurately characterise the impact of chaos in the Aquar¬ 
ius haloes themselves, but to obtain reasonable descriptions 
of these numerically simulated DM haloes to reflect in our 
results the expected stochasticity due to the morphological 
properties of this galactic component. 

As shown by Valluri et al. (2012]), taking initial con¬ 
ditions from a self-consistent model and evolving them in 
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Table 1. Main properties of the five Aquarius haloes at z = 0 from Springel et al.| ( |2008b[ ). The first column 
labels the simulation. From left to right, the columns give the virial radius of the DM halo, rooo; the virial 
mass, M200; the number of particles within 7*200, Abo a ; the particle mass, m p ; the concentration parameter, 
Cnfw; the intermediate to major, b/a and the minor to major, c/a, principal axial ratios computed using 
DM particles located within 6 to 12 kpc; the total stellar halo mass, M* (our stellar halo masses also 
includes the mass assigned to the bulge component in|Cooper et al.f2010) and the half-light radius from 


Cooper et al. 
Masses are in 


12010 ), 

Hgpdistances in kpc and velocities in km s —1 


ri/ 2 - 


Name 

^200 

M200 

[10 12 ] 

A200 

[ 10 6 ] 

m p 

[10 3 ] 

cnfw 

b/a 

c/a 

[10 s ] 

ri/2 

Aq-A2 

245.88 

1.842 

135 

13.7 

16.19 

0.65 

0.53 

3.8 

20 

Aq-B2 

187.7 

0.8194 

127 

6.4 

9.72 

0.46 

0.39 

5.6 

2.3 

Aq—C2 

242.82 

1.774 

127 

14.0 

15.21 

0.55 

0.46 

3.9 

53 

Aq-D2 

242.85 

1.774 

127 

14.0 

9.37 

0.67 

0.58 

11.1 

26 

Aq-E2 

212.28 

1.185 

124 

9.6 

8.26 

0.67 

0.46 

18.5 

1.0 


a slightly different potential should increase the fraction of 
chaotic orbits in the sample. Thus, our results are likely to 
overestimate the role of chaotic mixing in the systems under 
study. 


In |Navarro et al.j (1996 1997|, the authors introduced a 
spherical density profile that provides a reasonable fit to the 
mass distribution of DM haloes of galaxies in a very wide 
range of mass and redshift. However, it is now known that 
in a A-CDM cosmology DM haloes are not spherical as as¬ 
sumed by this potential. Instead, these are expected to be 
strongly triaxial and, furthermore, their shape is expected 


to vary as a function of galactocentric distance (see e.g. All 


good et al. 2006 Vera-Ciro et al. 20111. Introducing into 


our analysis the triaxiality of the galactic potential is of key 
importance, as this is one of the main sources of chaotic be- 


haviour together with the cuspy profile 1 

Voglis et al. 

2002 

Kandrup & Siopis|20031 Kalapotharakos et al. 2004; IV 

uzzio 

et al. 2005). 

triaxial extension 
axiality and radial 

of t 

Vogelsberger et al. (2008) presented 

lis profile that takes into account tri 


variation in shape. The associated potential, <I>tri, can be 
described by 


* raj = -£l n ( 1 + ^) , ( 1 ) 

where A is a constant defined as 

A _ _ G M 200 _ 

In (1 + cnfw) — cnfw/ (1 + cnfw) ’ 

with G the gravitational constant, M 200 the virial mass of 
the DM halo and cnfw the concentration parameter; r s = 
D20o/cnfw is a scale radius with 7*200 the virial radius. The 
triaxiality of this potential is introduced through r p , 

(r 3 + r)r e 

r p = ---, 

r a + r e 

where r is the usual galactocentric distance and r e an ellip¬ 
soidal radius defined as 


ratios and directions of the principal axes were computed 
using DM particles located within 6 to 12 kpc. Their values 
are listed in Table [T] 

Note that, in this approximation used to represent the 
underlying potential of DM haloes, the potential shape 
changes from ellipsoidal to near spherical at the scale radius, 


r a . Thus, for r -C r a , r p r e and for r > r s , r p ^ r (Vogels- 


|berger et al.|2008| ). We find that, to the 10% level, these ap¬ 
proximated analytic potentials can reproduce the true grav¬ 
itational potentials within the relevant distance range, i.e., 
r < 100 kpc. This is in agreement with the recent results 
presented by Bonaca et al. (20141, using the Via Lactea II 


simulation ( Diemand et al.|2008 l. 

The potential $tri admits, for r p < r a the power series 
expansion 


$ TM 


(- 1 ) 


n+1 



( 2 ) 


so it is analytic everywhere, and the condition r p < r a im¬ 
plies that r, r e < r a . 

Under the above assumption, r p /r a could be approxi¬ 
mated, up to r 2 /r 2 , by 


1 + 


1 - 


In spherical coordinates (r, 1 ?, 95 ), introducing the parameters 


£1 


^-1 

b 2 


£2 


C-1 


retaining terms up to r p /r a in | 2 | and neglecting a constant 
term, the potential takes the form 


4?tri( r*, $, « 4>o(r) + <&i(r) {(e 2 - ei)cos2$- 

— £1 cos 2<p + cos 2(i} + <fi) + cos 2($ — <p)}, (3) 

where 


W(i)’+«)■+©■• 

The quantities b/a and c/a represent the intermediate to 
major and the minor to major principal axial ratios and are 
defined such that a 2 + b 2 + c 2 = 3. In all simulations, the 


‘ho O') = 

Ar 

2ar a 

K 

'l’i (?•) = 

Ar 2 

2 arf 

K 


1 -— ) + (£1 + £2)^1 ( 0 , 

ar s 


2 

ar s 





















































On the relevance of chaos 5 


This approximation will be used in Section |4.2| when 
discussing diffusion. 


initial time to and a stopping time tf, we can define the 
OFLI as: 


2.3 Cosmological motivated initial conditions 


To investigate the efficiency of chaotic mixing on halo stars 
in the vicinity of the Sun we first need to model their distri¬ 
bution in phase space. Rather than stochastically sampling 
the phase space distribution associated with the potential 
presented in Section 2.2 Eq. |l]), we select from each halo 
the stellar particles within spheres centred at 8 kpc from 
the corresponding galactic centre. Following |Gomez et al.| 
1 2010), we choose for the spheres a radius of 2.5 kpc. This ra¬ 


dius approximately corresponds to the distance within which 
the astrometric satellite Gaia will be able to measure with 
high accuracy positions and velocities of an extremely large 
number of stars. As the final configuration of the five host 
DM haloes is strongly triaxial, we have rotated each halo to 
its set of principal axes and placed the corresponding local 
spheres along the direction of the major axis. This allows 
a direct comparison between the different haloes. As shown 
by G13, varying azimuthally the location of our spheres re¬ 
sults in local stellar densities that are, in general, an order 
of magnitude smaller than the observed value in the Solar 
Neighbourhood. Furthermore, the differences in number of 
resolved stellar streams within spheres located at different 
azimuthal angles mainly reflect changes in the local stel¬ 
lar density. Thus, we do not expect our particular choice of 
location for these spheres to affect significantly our results 
concerning the dynamical nature of Solar Neighbourhood¬ 
like volumes. Finally, we will only consider stellar particles 
that originally were members of accreted satellite galaxies. 
Any stellar particle associated with in-situ star formation 
are disregarded. 


2.4 Chaos Indicator: The Orthogonal Fast 
Lyapunov Indicator 

Now that the model and the volumes of interest have been 
introduced, the main goal of this section is to present briefly 
the preferred Cl used in the analysis: the OFLI (jFbuchard 
et al.j 2002 I, a particular variant of the FlQ 

Given an A r -dimensional Hamiltonian H. If we follow 
the time evolution of a unit deviation vector w (t) for a given 
solution of the equations of motion 7 (t), initially chosen nor¬ 
mal to the energy surface (i.e. in the direction of VH, see 
Barrio| 2015 I, take its orthogonal component to the flow at 
time t, w(t) x € R, and retain the largest value between an 


2 Even though we adopt the OFLI as our primary chaos indicator, 
this study is supported by similar results based on other CIs, 
such as the LI (a short introduction to the basic idea behind 
CIs and a definition of the most popular Cl, the LI, are given 


in Appendix [A} , the MEGNO 

0incotta & Simo 2000j Cincotta 

et al.||2003| |Gozdziewski et al. 1 

2005| |Compere et al.||2011 

), the 

GALI ('Skokos et al. 2007 200£ 

Manos & Athanassoula 

2011 

Manos et al.|2012b and the RLI 

Sandor et al.| 2000) )2004 

2007 


Szell et al. 2004). Thus, the orbital classification obtained with 
the approximate galactic potential described in Section [2.2| is very 
robust. However, for the sake of brevity the results based on these 
other CIs are not presented in this work. 


0FLI 7 (t/) = sup 


t 0 <t<tj 

for the orbit 7 . The OFLI 7 tends to infinity as time increases 
for both non-periodic regular and chaotic orbits. The growth 
of OFLI 7 is exponential with time if 7 is a chaotic orbit. 
The OFLI grows linearly with time for resonant and non¬ 
resonant regular orbits (on a logarithmic scale), with differ¬ 
ent rates, and oscillates around a constant value for periodic 


orbits (for further details we refer the reader to Froeschle et 
ak||1997| |Froeschle fc Legal |2000[|Legafc Froeschle||2001| 


Fouchard et al.|2002 Guzzo et al.|2002 Barrio||2015 \. 

In what follows, we integrate the orbits and compute 
the CIs using the LP-VIcode program. LP-VIcode is a fully 
operational code which efficiently calculates a suite of 10 CIs 
in any number of dimensions (see Carpintero et al. 2014 1 . 
The hardware we used for these experiments was an Intel 
Core i5 with four cores, CPU at 2.67 GHz and 3 GB of 
RAM. The version of the gcc gfortran compiler was 4.8.2. 


3 BUILDING UP A BASIC UNDERSTANDING 


Chaos indicators are used by dynamicists to identify and 
characterise the interplay between regular and chaotic com¬ 
ponents of diverse dynamical systems. For instance, CIs are 
a popular means of quantifying the impact of chaos on the 
dynamical evolution of self-consistent stellar systems, us¬ 
ing non-evolving analytic models of the underlying poten¬ 


tial (Kalapotharakos & Voglis 2005. Manos & Athanassoula 


2011 Zorzi fc Muzzio||2012 1 . These CIs, such as the OFLI, 

provide a reliable and straightforward way to estimate the 
chaos onset times of orbital sets, i.e. the timescale at which 
stars commonly associated with chaotic orbits will effec¬ 
tively reveal their chaotic behavioui^] I 11 the following ex¬ 
periment we apply the OFLI to a regular orbit, a sticky 
orbit and a chaotic orbit in order to show what should be 
expected from the indicator in each case. 

Sticky orbits are, for the purpose of this paper, chaotic 
orbits that behave regularly on timescales comparable to a 
Hubble time. In what follows, we classify an orbit as sticky if 
it exhibits chaotic behaviour only on timescales larger than 
10 Gyr. For a more formal definition of sticky orbits, please 
see Appendix [B| 

In Fig. |T]we present typical examples of regular orbits, 
sticky orbits and chaotic orbits in a triaxial extension of the 
NFW potential (Eq. ,0) with the parameters of the DM 
halo Aq-A2 (Table [I]). The (rounded) initial conditions for 
these orbits are presented in Table [2] We integrate simul¬ 
taneously the equations of motion (Eqs. (Al I) and the first 


variational equations (Eqs. (A2 1 ). The regular and sticky or¬ 
bits are integrated for a rather long integration time of 1000 
Gyr, while the chaotic orbit is integrated over only 10 Gyr, as 
this timescale is enough to characterise its behaviour. The 
time-step used for the regular and sticky orbits is 1 Myr, 
and 0.01 Myr for the chaotic orbit. The (rounded) binding 


3 Note that the chaos onset time is not the same as the commonly 
used Lyapunov time. 
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Regular orbit: 950-1000 [Gyr] 


Sticky orbit: 890-1000 [Gyr] 


Chaotic orbit: 5-10 [Gyr] 



Figure 1. Examples of regular orbits (left panels), sticky orbits (middle panels) and chaotic orbits (right panels) in the triaxial extension 
of the NFW model, for different time intervals. The fact that different volumes are occupied by the orbits in different time intervals (see 
the middle and right panels) is an indication of chaotic behaviour. 


Table 2. Initial conditions and binding energies for our examples of regular, sticky and chaotic orbits. 
Distances are in kpc, velocities in km s 1 and binding energies in km 2 s —2 . 


Type of orbit 

X 

y 

z 

v x 

Vy 

V z 

E 

Regular 

8.219 

-0.652 

-2.203 

-5.795 x 10“ 3 

102.95 

-4.745 

-217691.706 

Sticky 

5.865 

0.263 

-0.346 

239.350 

333.547 

57.208 

-152469.329 

Chaotic 

5.731 

0.531 

0.443 

6.029 

-0.366 

23.691 

-238489.404 


energies (E) of the three orbits are included in Table [5] The 
integrator conserves energy to an accuracy of one part in 
10“ 12 or less for all the experiments throughout the paper. 

In Fig. [T] we show only the first and last intervals of 50 
and 110 Gyr for the regular and the sticky orbits, respec¬ 
tively, to illustrate the different behaviours that characterise 
these types of orbit in configuration space. The behaviour 
of the chaotic orbit is presented for two consecutive 5 Gyr 
intervals. 

From Fig. [l] it is clear that the regular orbit has a very 
similar shape in the two different time intervals, even though 
these are separated by 900 Gyr. Much more significant dif¬ 
ferences are apparent in the trajectories of the sticky orbit 
sampled at widely separated time intervals (0 to 110 Gyr for 
the middle top panel and 890 to 1000 Gyr for the middle bot¬ 
tom panel). This illustrates the chaotic nature of the sticky 
orbit. Nevertheless, it takes the orbit more than 8 Hubble 
times (roughly 110 Gyr) to exhibit its true nature. In the 
rightmost panels we show a chaotic orbit for two consecu¬ 
tive 5 Gyr intervals. Even in these much shorter intervals, 
the evolution in the trajectory is evident. 

In Fig. [2] we present the characteristic behaviour of the 
OFLI for the three types of orbit showed in Fig. |T] In this 
case, all three orbits have been integrated using a time-step 
of 1 Myr. Notice the linear evolution of the indicator for 


regular orbits and the exponential growth corresponding to 
sticky and chaotic orbits. As expected from our previous 
discussion, for the sticky orbit, it takes the OFLI ~ 100 Gyr 
to start growing exponentially, time at which the chaotic 
behaviour of this particle is revealed. Instead, for the chaotic 
orbit, it only takes the OFLI a few Gyr to start showing an 
exponential growth. 


4 THE ACTUAL RELEVANCE OF CHAOS: 
SOLAR NEIGHBOURHOOD-LIKE 
VOLUMES 


In order to characterise the impact of chaotic mixing on 
the phase space distribution of the Solar Neighbourhood, in 
this section we examine two central points: i) the distribu¬ 
tion of chaos onset times for stellar particles within Solar 
Neighbourhood-like volumes and ii) the rate of diffusion 
due to chaotic mixing, a mechanism that can lead to large 
variations of the integrals of motion. Our goal is to explore 
whether chaotic mixing can be strong enough to erase sig¬ 
natures of merger events in the neighbourhood of the Sun. 

To tackle (i) we select particles from the five Aquar¬ 
ius stellar haloes in Solar Neighbourhood-like volumes (see 
Section 2.31. We then quantify the fraction of particles on 
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Figure 2. Time evolution of the OFLI for the three orbits intro¬ 
duced in Fig. [T] Notice the logarithmic scale on both axes. The 
exponential growth of the indicator for chaotic motion is clearly 
observed for the chaotic and sticky orbits. 


regular, sticky and chaotic orbits. We do this by means of 
the OFLI, which allows us to estimate efficiently the distri¬ 
bution of chaos onset times and so identify the characteristic 
timescales over which chaotic mixing becomes relevant. To 
address (ii), we measure the diffusion of pseudo-integrals of 
motion for large ensembles of test particles that are initially 
nearby in phase space. 


4.1 The importance of timescales 

First we look for a relationship between the expected OFLI 
behaviours of different orbit types (described above) and 
more readily interpreted measures of the local (stream) den¬ 
sity. After this, we estimate chaos onset times to characterise 
the actual amount of chaos manifested on physically mean¬ 
ingful timescales. 


4-1.1 The time evolution of the OFLI for initially nearby 
particles in phase space 

In the previous section we have shown that our triaxial DM 
haloes admit a wide range of different orbital behaviours. In 
particular, we have shown an example of a very chaotic orbit 
with a chaos onset time shorter than 5 Gyr. These orbits can 
potentially play an important role in shaping the present- 
day phase space distribution of our Solar Neighbourhood. 
As shown by Vogclsbcr ger et al.| ( |2008| ), the local density 
in the neighbourhood of a particle moving on a chaotic or¬ 
bit decreases exponentially with time. As a result, stellar 
streams moving on such orbits will experience a rapid phase 
space mixing process which can erase signatures of their ac¬ 
cretion history. In contrast, the local density of a particle 
moving on a regular orbit decreases as a power law function 
of time, with exponent less than or equal to 3 (a signifi- 
cantly lower rate, |Helmi fe White 1999. Vogclsbcrger et ah] 


|2008||G6mez et al.||20l 0 l Hence, the probability of finding 

streams on regular orbits in the Solar Neighbourhood may 
be somewhat higher. 

With these ideas in mind, G13 followed the time evolu¬ 
tion of the local density of accreted stellar particles that, at 


2 = 0, were located within different Solar Neighbourhood¬ 
like volumes. The goal was to explore whether particles that 
appeared to be smoothly distributed in phase space, and 
thus not associated with any resolved stellar stream, were 
on chaotic orbits. If not, the lack of clustering in phase space 
for these particles could be due to the limited numerical reso¬ 
lution of the simulations. G13 fit a power law function to the 
time evolution of the local density around each star particle 
and determined the rate at which this local density decreases 
with time. They assumed that a power law fit to the local 
density of stream particles on a chaotic orbit should yield 
an exponent greater than 3. Unfortunately, due to the finite 
resolution of their simulations, the local volumes used to 
track the time evolution of density in G13 were rather large 
- namely, spheres of radius equal to half of the apocentre of 
the particle’s orbit. Such large volumes are not problematic 
if the goal is to detect resolved stellar streams. However, it 
is not clear that such fits reflect the true dynamical nature 
of the underlying local stream densities. 

In the following experiments we explore the connection 
between the OFLI and the time evolution of the local density 
of a stream. We also explore the effects of the size of the local 
volume on the characterisation of the dynamical nature of a 
stream through the time evolution of its local density. 

Figure [3] shows the trajectories of two specific stellar 
particles in the Aq-A2 DM halo. These particles were cho¬ 
sen as clear examples of sticky and chaotic orbits, and at 
2 = 0 are found in a sphere of radius 2.5 kpc located at 
8 kpc from the centre of the halo at 2 = 0. We refer to 
these stellar particles as guiding particles A and B; their es¬ 
timated chaos onset times are ~ 113 and 5 Gyr, respectively. 
To relate the OFLI to the time evolution of local density 
of streams, we distribute an ensemble of 1000 massless test 
particles in a small phase space volume around both guiding 
particles. These volumes are defined by multivariate Gaus- 
sians in phase space with initial dispersions cr x = 0.2 pc and 
cr v = 1 km/s. Test particles are initially distributed such 
that their maximum separation with respect to the guiding 
particle is less than or equal to \/3cr x and \/3 <t v respectively. 
As we show below, such small phase space volumes are nec¬ 
essary to characterise accurately the dynamical behaviour 
of the local density around our particles A and B. 

The ensembles of test particles and their corresponding 
guiding particles are then integrated in our triaxial NFW 
model (Aq-A2 DM halo parameters as Table [T]) for 10 Gyr, 
with a timestep of 1 Myr. During the integration, the OFLI 
is computed for each particle (test and guiding). 

The left and middle panels of Fig. [4] show the time evo¬ 
lution of the OFLI for the guiding particle A and its asso¬ 
ciated ensemble of test particles, eA. Guiding particles are 
depicted in cyan and test particles in black. The time inter¬ 
vals shown in these panels are the same as those in Fig. [3] i.e. 
0 to 10 Gyr (left panel) and 110 to 120 Gyr (middle panel). 
Note that separate ensembles of test particles are sampled 
from the Gaussian kernel defined above at the start of each 
interval shown. 

On the left panel of Fig. [4] we show that the OFLI for 
particle A increases linearly with time over the first 10 Gyr 
of evolution, indicative of a regular orbit as described in 
Section [3] The indicator shows similar behaviour for all the 
test particles eA over the same interval (black solid curves). 
This indicates that the initial distribution of test particles 
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Particle A: 0-10 [Gyr] 


Particle A: 110-120 [Gyr] 


Particle B: 0-10 [Gyr] 



Figure 3. Trajectories for both guiding particles (A and B) over different time intervals. The leftmost panel shows the trajectory of 
particle A from 0 to 10 Gyr and the central panel its trajectory from 110 to 120 Gyr. The rightmost panel shows the trajectory of particle 
B from 0 to 10 Gyr. The orbit associated with particle A has a chaos onset time larger than 110 Gyr. Thus, its shape only starts to 
change after that time. In case of particle B, the chaos onset time is shorter than 10 Gyr. 
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Figure 4. Time evolution of the OFLI for both guiding particles (A and B) and their corresponding test particle ensembles (eA and 
eB). Guiding particles are shown in cyan and test particles in black. Left panel: OFLI for particle A and ensemble eA in the interval 
0 to 10 Gyr. Middle panel: OFLI for particle A and ensemble eA from 110 to 120 Gyr. Right panel: OFLI for particle B and ensemble 
eB in the interval 0 to 10 Gyr. Notice the logarithmic scale. The left and middle panels show the typical behaviour of the phase space 
surrounding a sticky orbit: regular behaviour followed by exponential growth of the OFLI at the chaos onset time. The last panel shows 
a much more rapid onset of exponential growth due to the high degree of chaos in the phase space neighbourhood of particle B. 


accurately samples the phase space volume of the guiding 
particle. 

The middle panel of Fig. [4] shows results for another 
time interval, from 110 to 120 Gyr, for the same guiding 
particle A. The exponential growth of the OFLI in this in¬ 
terval (starting at ~ 113 Gyr) implies a chaotic orbit. As 
before, we find similar behaviour of the OFLI for all the test 
particles in the corresponding ensemble, including the onset 
time of exponential growth. Clearly, particle A moves on a 
sticky orbit. 

The right panel of Figure [4] shows the time evolution 
of the OFLI for another guiding particle, B, with a chaotic 
orbit. Note that the chaotic nature of this particle’s orbit 
becomes apparent on a much shorter timescale, ~ 5 Gyr. 
As in the second time interval shown for particle A, the 
vast majority of test particles in ensemble eB show rapid 
growth of the OFLI, reflecting the behaviour of their guiding 
particle. However, we find a much larger spread in the rate 
of the exponential growth in the case of particle B than in 
the case of the sticky orbit of particle A. 

We have demonstrated that test particles distributed 
within initially small phase space volumes around a guiding 
particle show very similar behaviour of the OFLI. This indi¬ 
cates that the time evolution of an ensemble of such particles 
can be used as an indirect indicator of the nature of any guid¬ 


ing particle. As previously described, this idea was exploited 
by G13 as a means of quantifying the nature of stellar orbits 
within Solar Neighbourhood-like volumes in the Aquarius 
simulations. However, the local stream density around each 
stellar particle in G13 was computed by following the neigh¬ 
bouring stellar particles that, at the corresponding stream’s 
formation time, were located within spheres of radii larger 
than or equal to 4 kpc radius, much larger than the kernel 
size we use in Figure [4] We now explore whether the above 
result holds even when larger local volumes are considered 
for the initial distribution of test particles. 

For this experiment we consider a guiding stellar parti¬ 
cle that is moving on a regular orbit. In Figure [5] we show 
the time evolution of the OFLI for test particles initially dis¬ 
tributed over configuration space kernels of different sizes. In 
all cases the particles were distributed over the same Gaus¬ 
sian kernel in velocity space, with a v = 1 km/s. The left 
panel shows the results obtained for <r x = 0.005 kpc. Within 
this relatively small sphere, the evolution of the OFLI for 
most test particles still reflects the behaviour of the guiding 
particle. 

In the middle panel of Figure [5] we use a configura¬ 
tion space kernel with cr x = 0.1 kpc. Although the majority 
of test particles still show a regular behaviour in this case, 
we start to find a significant fraction that show a chaotic 
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Figure 5. Time evolution of the OFLI for 1000 test particle initial conditions sampled around a guiding particle on a regular orbit. The 
radius of the initial sphere is indicated. Notice the logarithmic scale. The bigger the sphere around the guiding particle, the more diverse 
the behaviours of the neighbouring particles. 


behaviour in the same integration period. The right-hand 
panel shows results for cr x = 4 kpc. It is clear from this 
panel that such an extended initial distribution of test par¬ 
ticles does not accurately reflect the dynamical behaviour 
of the guiding particle. The results presented in G13 may 
overestimate the fraction of chaotic orbits within each Solar 
Neighbourhood-like volume. 


4-1.2 Connecting the OFLI to a measure of the local 
stream density 

We have shown that initially nearby particles in phase space 
show similar time evolution of the OFLI, which implies a 
common dynamical behaviour. We will now explore whether 
the behaviour of the OFLI accurately reflects the time evolu¬ 
tion of local density along a stream. For this purpose, we fol¬ 
low the evolution of the local density around a sticky stellar 
particle and a chaotic stellar particle for periods of 10 Gyr. 
We refer to these stellar particles as guiding particles. Both 
guiding particles are located in a Solar Neighbourhood like 
sphere at 2 = 0. Estimates of the chaos onset times for these 
guiding particles are ~ 80 and 3 Gyr, respectively. Note that 
the following results do not depend strongly on our specific 
choice of stellar particles. 

As before, we distribute ensembles of test particles 
around each guiding particle. The test particles are initially 
distributed as explained in Section Em with <T X = 0.2 pc 
and ov = 1 km/s. However, to accurately track the time 
evolution of the local density around both guiding particles 
for periods of 10 Gyr, a larger number of test particles, 10 4 , 
is considered for each ensemble. Two different time intervals, 
separated by several Hubble times (~ 113 Gyr), are consid¬ 
ered. The integration timestep is 0.1 Myr. As expected from 
their chaos onset times, during the first time interval the 
sticky orbit behaves like a regular orbit while the chaotic 
orbit shows its true nature. During the second time interval 
the sticky orbit behaves like a chaotic orbit. To estimate the 
local density at every timestep, we count the number of test 
particles within a radius of 0.1 kpc around the guiding parti¬ 
cle, and also discard from further consideration any particles 
beyond a radius of 2 kpc. The number of test particles within 
the 0.1 kpc sphere is then normalised by the initial number 
of test particles. We call this quantity the normalised num¬ 
ber of neighbouring particles: N;. Here i = S, C refers to the 


normalised densities associated with the sticky and chaotic 
orbits, respectively. 

On the top panel of Fig. [6] we present, with dark grey 
dashed and green solid lines, the time evolution of Ns for 
both time intervals (labelling the first 10 Gyr as interval 1 
and the second 10 Gyr time interval as number 2). On the 
same panel we also present, with a red solid line, the time 
evolution of Nc for the first 10 Gyr time interval. It is clear 
from this panel that the local density around the sticky orbit 
decreases in time significantly more slowly when the guid¬ 
ing particle’s orbit is approximately regular. Notice that, in 
this regime, we are able to track the evolution of the local 
density for the full 10 Gyr period of integration. Conversely, 
in the chaotic regime, the number of neighbouring parti¬ 
cles within 0.1 kpc, that have never been further than 2 
kpc, becomes insufficient to track the local density after ~ 6 
Gyr. Notice that the sticky particle and the chaotic particle 
show a similar evolution of their local densities during the 
chaotic regime. The more rapid decay of Nc reflects the more 
chaotic nature of this orbit. On the bottom panel of Fig. [6] 
we show that a power law function describes well the time 
evolution of the local density around the sticky particle dur¬ 
ing the first 10 Gyr time interval; specifically Ns oc t -1 ' 3 , 
Conversely, an exponential relation is required to describe 
the time evolution of the local density around the chaotic 
particle Nc oc e -1 ' 46t . 

This analysis demonstrates a very strong connection be¬ 
tween the evolution of local (stream) density around a given 
guiding particle and the characterisation of its orbit pro¬ 
vided by the OFLI. If the OFLI shows a linear growth, indi¬ 
cating regular behaviour, then the local density is expected 
to decrease as a power law with index less than or equal to 
3. Likewise, an exponential growth of the OFLI reflects an 
exponential decay of the corresponding local density. 


4-1.3 The distribution of chaos onset times within Solar 
Neighbourhood-like volumes 

In the previous section we have shown that the OFLI is a 
powerful tool to characterise the time evolution of the local 
(stream) density around any stellar particle in our simula¬ 
tions. We will now use this indicator to quantify robustly the 
fraction of stellar particles within Solar Neighbourhood like 
volumes that are moving on regular, sticky and chaotic or¬ 
bits. If chaotic orbits are common in a phase space volume 
















10 


N. P. Maffione et al. 



0 2 4 6 8 10 


Time [Gyr] 



0 2 4 6 8 10 

Time [Gyr] 


Figure 6. Top panel: time evolution of the normalised number 
of neighbouring particles for the sticky orbit, Ng, over two non- 
consecutive 10 Gyr time intervals. The first interval (labelled 1) 
is taken when the guiding particle moves on a regular orbit and 
the second interval (labelled 2) is taken when the guiding particle 
moves on a chaotic orbit. We also show the time evolution of 
the normalised number of neighbouring particles for the chaotic 
orbit, Nq, over the first time interval. Bottom panel: a power law 
fit for Ng and an exponential fit for Nq for the first interval only. 
Notice the logarithmic scale. The local density of stellar streams 
decreases with a power law function along a regular orbit and at 
an exponential rate along a chaotic orbit. 


local to the Sun, then many phase space substructures (aris¬ 
ing for example from past accretion events) may have been 
erased, due to the much shorter mixing timescales associated 
with such orbits. 

In this experiment we will consider five different Solar 
Neighbourhood-like volumes, each extracted from a differ¬ 
ent Aquarius stellar halo. In all cases the spheres are cen¬ 
tred at 8 kpc from the galactic centre and have a radius 
of 2.5 kpc (see Section [2.3| ). In order to compute the OFLI 
for each stellar particle within these spheres, we integrate 
the equations of motion (Eqs. (Al l) together with the first 
variational equations (Eqs. (A2|), assuming smooth triax- 
ial NFW DM haloes (Eq. (JX I) with parameters as given in 
Table [I] The timestep of integration is 1 Myr and the to¬ 
tal integration time 1000 Gyr. We use such a long timespan 
(more than 100 times the likely age of the Milky Way) to 
identify very sticky orbits reliably. Finally, we compute the 


distribution of estimated chaos onset times to obtain the 
percentages of particles moving on sticky and chaotic orbits. 

As previously mentioned in Section [3] we consider as a 
sticky orbit any orbit that has a chaos onset time larger than 
10 Gyr (a Hubble time, roughly speaking). This definition 
is arbitrary. However, it allow us to make a clear distinction 
between orbits that could show some degree of chaotic mix¬ 
ing within a physically meaningful timescale and those for 
which chaos is completely irrelevant. 

As described in Section |2.3| only accreted stellar par¬ 
ticles are analysed in these experiments. Table [3] lists the 
total number of stellar particles considered in each stellar 
halo, N*. This number is slightly different from the number 
of stellar particles quoted in G13 (see their Table 3). This is 
due to a small number of stellar particles, N°, that needed 
some numerical readjustments during the integration - the 
computation of these particles was stopped in order to keep 
the computing time of the integrator bounded. 

The orbits of the particles selected for further study 
were then classified according to the shape of of their OFLI 
time evolution curve. Individually inspecting each curve to 
estimate the chaos onset time is unfeasible for large samples 
of orbits. We therefore introduce a threshold OFLI value 
based on an upper limit for the typical linear behaviour seen 
for regular orbits. This threshold evolves linearly with time 
and envelopes all the curves that present a linear behaviour. 
On Figure [7] we show an example of this procedure. The dif¬ 
ferent lines show the time evolution of the OFLI for 1388 
particles located within the Aq-A2 Solar Neighbourhood¬ 
like sphere. The threshold is indicated with a blue solid line. 
Every time the OFLI of any given stellar particle crosses 
this threshold, the corresponding particle is classified as ei¬ 
ther sticky or chaotic. We reject all threshold crossings in 
the first Gyr of evolution, because this period can be con¬ 
sidered as a transient stage of the indicator. Chaotic orbits 
are defined by threshold crossing within the first 10 Gyr of 
their evolution. This 10 Gyr barrier is depicted in Fig. 0by 
a vertical dashed blue line. We denote the threshold crossing 
times of each stellar particle by T* (our estimation of the 
chaos onset time). The different colour-coded curves in Fig¬ 
ure [7] show the results of this classification applied to halo 
Aq-A2. Similar results are found for the other four haloes, 
which we do not show for the sake of brevity. Note that, in 
a very small number of cases, the OFLI of an orbit crosses 
the threshold early on, but later continues to evolve linearly 
with time. As a consequence, the number of chaotic orbits 
found within each volume may be slightly overestimated. 

In Table [3] we summarise the result of this experiment. 
It is interesting to see that, even in significantly triaxial and 
cuspy potentials, a significant fraction of orbits are regular 
even after 1000 Gyr of evolution. In all haloes, we find that 
« 30% of the stellar particles are moving on regular orbits. 
It also striking to find that, in all haloes, the fraction of 
stellar particles on sticky orbits is larger than 45%. As a 
consequence we find that only < 20% of orbits could be 
experiencing some degree of chaotic mixing, regardless of 
the scale and shape of the halo. 

Stellar particles living in the innermost regions of each 
stellar halo are likely to be more bound and to have shorter 
dynamical timescales than those populating the outer galac¬ 
tic regions. It is interesting to study whether the distribu¬ 
tion of orbital dynamical timescales (hereinafter Dyt s ) plays 
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Table 3. Dynamical distributions of particles at z = 0 from the Aquarius Project. The first column labels 
the simulated DM haloes. From left to right, the columns give: the number of particles, N*; the number of 
particles which computation was stopped, N °; the number of particles on regular orbits, TV^; the number of 
particles on sticky orbits, the number of particles on chaotic orbits, Nq] the number of particles with 
values of the Dyt s < Dy™ ax , N£ and the medians of Dyt s for the sticky, /xjy 2 , anc ^ the chaotic distributions, 

m f /2 in Gyr. 


Name 

N* 

N° 


N s * 

N* 

N* 

,,S 

"l/2 

1J C 

Pl/2 

Aq-A2 

1388 

12 

437 (31.48%) 

644 (46.4%) 

307 (22.12%) 

1310 (94.38%) 

0.295 

0.256 

Aq-B2 

10385 

355 

3618 (34.84%) 

5069 (48.81%) 

1698 (16.35%) 

10273 (98.92%) 

0.417 

0.396 

Aq—C2 

2291 

43 

756 (31.62%) 

1015 (46.63%) 

520 (21.75%) 

2275 (99.3%) 

0.297 

0.263 

Aq-D2 

3745 

108 

1087 (29.03%) 

2213 (59.09%) 

445 (11.88%) 

3679 (98.24%) 

0.518 

0.409 

Aq-E2 

1877 

80 

589 (31.38%) 

888 (47.31%) 

400 (21.31%) 

1853 (98.72%) 

0.494 

0.426 



Time [Gyr] 


Figure 7. Time evolution of the OFLI for the 1388 particles 
considered in the Aquarius Project for the Aq—A2 DM halo and 
within an interval of time long enough to identify very sticky 
orbits (1000 Gyr). The upper limit used as a threshold for regular 
motion is depicted in solid blue. The 10 Gyr threshold is depicted 
with a vertical dashed—blue line. Notice the logarithmic scale. 
The three orbital components, i.e. the sticky, the regular and the 
chaotic components, are clearly distinguished by using the OFLI 
with both simple thresholds. 

an important role in separating sticky from chaotic orbits. 
If chaotic orbits are preferentially found in the innermost 
galactic regions, then we may have been able to identify 
them simply because their corresponding Dyt s are small 
enough to reveal their nature in short integration times. On 
the other hand, for orbits probing the outer regions of the 
halo, with larger values of Dyt s , very long integration times 
could be required for an accurate characterisation. It is thus 
important to understand whether or not our quantification 
of regular, sticky and chaotic orbits is biased by differences 
in the relative distributions of their orbital timescales. 

To explore this we define a characteristic Dyt s as the 
time for a particle on a given orbit to undergo two changes 
of the sign of its velocity component along the major axis 
of the stellar halo. In Figure [8] we present, for halo Aq- 
A2, the distribution of Dyt s as a function of the estimated 
chaos onset time, T*. Regular, sticky and chaotic orbits are 
depicted in black, green and red, respectively. As before, 
this classification is based on their values of T*. It is very 
clear from this panel that chaotic and sticky orbits are not 


Figure 8. The threshold crossing time, T*, as a function of the 
dynamical timescale, Dy ts- The chaotic orbits are depicted in 
red, the sticky orbits in green and the regular orbits in black as 
throughout the paper. The vertical dashed blue line is the thresh¬ 
old for Dy™ ax , i.e. 0.631 Gyr for the halo Aq-A2. The distribu¬ 
tions for both chaotic and sticky orbits are very similar. 


strongly segregated in Dyt„. Note that many sticky orbits 
with T* ^ 100 Gyr can be found with values of Dyt s as 
small as ~ 0.2 Gyr. A similar result is observed for regular 
orbits. Even though sticky orbits show a tail towards large 
values of Dyt s , the populations of chaotic and sticky orbits 
seems to be characterised by similar distributions. This can 
be seen from the medians of the two distributions, given in 
Table 0 A vertical dashed-blue line indicates the location 
of the chaotic orbit with the highest value of Dyt s : Dy^ ax , 
which is 0.631 Gyr for halo Aq-A2. In Table [3] we present 
the fraction, N£, of orbits (of all classes) with values of Dyt s 
less than the corresponding Dy“ ax . In all haloes, N£ > 95%. 
In addition, we find that only ss 12% of the sticky orbits are 
above this threshold in halo Aq-A2. This number is reduced 
to < 2% for the other haloes. This means that short orbital 
timescales are not the dominant factor that distinguishes be¬ 
tween sticky and chaotic motion for stellar halo particles in 
our Solar Neighbourhood-like volumes. Instead, the results 
of this simple analysis demonstrate that the key distinction 
is the dynamical properties of the surrounding phase space 
volume. 

So far, we have shown that a small but non-negligible 
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fraction of orbits in Solar Neighbourhood-like volumes could 
indeed exhibit chaotic mixing. In what follows we will discuss 
the extent to which such mixing can erase nearby signatures 
of early stellar accretion onto the galaxy, within physically 
relevant periods of time. 


4.2 Global Dynamics and Diffusion 


4-2.1 Basic concepts 

In this section we discuss a mechanism that could lead to 
global chaotic mixing. In terms of the orbits of stars, roughly 
speaking, chaotic mixing means that all trajectories start¬ 
ing in a small neighbourhood of a given point in phase space 
(or on an energy surface) loose the memory of their initial 
conditions with time and eventually appear uncorrelated. 
For ^-dimensional systems (N ^ 3), since KAM tori no 
longer divide the energy surface (see Appendix [B] for a fur¬ 
ther description and references), it has been conjectured that 
any orbit lying in a thin chaotic layer around any resonance 


might visit the whole so-called Arnold web (Arnold 1964 
Chirikov| 1979). 


In his note, Arnold proved the existence of motion along 
the stochastic layer of a given resonance in a rigorous way, 
for a rather simple near-integrable Hamiltonian. He demon¬ 
strated that, for a very small perturbation it is possible to 
find a trajectory in the vicinity of the separatrix of a par¬ 
ticular resonance that connects two points separated by an 
arbitrarily large distance, i.e. independent of the size of the 
perturbation, on a very long timescale. Arnold’s proof rests 
on the existence of a chain of hyperbolic tori along this res¬ 
onance that may provide a path for the orbit - if these tori 
are very close to each other, an orbit could transit over that 
chain. Since every torus in the chain is labelled by an action 
value or unperturbed integral, a large but finite variation 
of this action could take place. This mechanism, which per¬ 
mits motion along the resonance stochastic layer, is known 
(in the mathematical literature) as the Arnold Mechanism, 
while the term Arnold diffusion generally refers (in the phys¬ 
ical literature) to a global phase space instability (for details 
Giorgilli|1990 Lochak|1999 Cincotta|2002| ), that is any 


(chaotic) orbit could visit the full Arnold web in a finite 
time. The problem of how to extend the Arnold mechanism 
to a generic Hamiltonian remains unsolved. One of the main 
difficulties is related to the construction of such a chain of 
tori. 

Regardless this severe limitation to understand Arnold 
diffusion as a global instability, it is assumed (in the physical 
literature) that Arnold diffusion does occur, and is responsi¬ 
ble for chaotic mixing (see for instance the discussion given 
in the last section of Cincotta et al. 2006). For instance, 


assuming that in the phase space of steady state galaxies 
the chaotic component is a well-connected region (through 
some type of diffusion), |Merritt| (j !999| l extended the clas¬ 
sical Jeans theorem, which was formulated only for regular 
non-resonant orbits ( Binney fe Tremaine|1987 1, to take into 
account chaotic motion. However, as we will show in this 
section, chaotic diffusion does not play any significant role 
in connecting the whole chaotic component. 

In spite of the mathematical difficulties in dealing with 
this conjecture as a global property, a local formulation 
shows that chaos needs to be considered in the limit when 


t —> oo in order to observe any significant variation of the un¬ 
perturbed integrals. This suggests tha t (strict) Arn old diffu- 
sion is irrelevant in the real worlcj^] |Chirikov & Vecheslavov 
1993 Ci ncotta|2002 I. 


In those real systems exhibiting a divided phase space, 
where the chaotic component is relevant (i.e. has a positive 
non-negligible measure) the timescale for any diffusion (not 
Arnold diffusion) would be much shorter but still very long 
(see for instance 
C incotta| 20041. 

In the particular model considered here, we claim 
that chaotic mixing is almost irrelevant on cosmological 
timescales. Although this is the aim of a forthcoming paper, 
we would like to show, by simple arguments and computa¬ 
tions, that the so-called chaotic diffusion does not work for 
the model here consider. 


Chirikov fe Vecheslavov| 1997 Giordano & 


4-2.2 Analytic description 


We have already shown that $tri can be approximated by 
Eq. (§ (see Section [2.2[), which for simplicity we recast as 


4>tri « 4> 0 (r) + $ 1 (r)[/i 1 cos2ip + p 2 cos2t?+ 

+ p -3 cos2($±<p)], 

where p a = e 2 — £1, —£1, £i/2, and since £1 and e 2 are as¬ 
sumed to be small parameters, /i s <l. 

Therefore, the Hamiltonian takes the form: 

?f(p,r) = H 0 (p,r,d) + 4>i(r), 


with 


'Wo(P ,r,9) 


and 


d p|_ P% 

2 2r 2 2r 2 sin 2 d 


+ $o (r), 


<f>i(r) = 4>i(r) [p\ cos2</> + p 2 cos2i? + pz cos 2(d ± ip )\, 
where 

p r = r. Pi} = r 2 d, p v = r 2 ip sin 2 d. 


In fact, Ho is an integrable Hamiltonian being 

Ho = E 0 , L z =Pcp, L 2 = pi + p 2 esc 2 d, 

the three global unperturbed integrals, while 4>i can be con¬ 
sidered as a small perturbation. The dependence of $1 on 
(d, p) leads to variations of the angular momentum and its 
components. Indeed, 


dL z 

dt 

dl? 

dt 




2p v dig 
sin 2 d dp 


which are of order p B and therefore assumed to be small. 

We are aware that the assumption /1 , ~£i,E 2 < 1 holds 
only marginally for any of the DM halo models considered 
here. However, the approach described above provides an 
appropriate physical insight into the problem. 


4 Further discussion about this instability and the connection 
between the mathematical and physical approach can be found 
in 1 


Guzzo & Lega (2013); 


Cincotta et al. 


(2014). 
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Figure 9. Ranges in L 2 (logarithmic scale) and L z for the whole 
set of 1400 particles (grey) of the Aq-A2 DM halo. In black, the 
region of the plane to be considered in the experiments. 


4-2.3 Numerical experiments 

In order to perform numerical experiments on a given con¬ 
stant energy surface on the (L 2 ,L 0 ) plane that mimics our 
Solar Neighbourhood, we consider the triaxial extension of 
the NFW model (Eq. 0 ) with parameters corresponding to 
the Aq-A2 DM halo. Initial conditions for an ensemble of 
test particles are obtained as follows. 

First, we fix ( xo,yo,zo) = (8,0,0) kpc (i.e. the position 
of the Sun), and adopt the mean value of the energy distribu¬ 
tion of the stellar particles located within a 2.5 kpc sphere, 
(E) = Eo ~ —204449 km 2 s -2 . The remaining two phase 
space coordinates are obtained by sampling a region of the 
(L 2 , L^) plane with a regular grid. The region of the (L 2 , L z ) 
plane explored is determined such that, in both dimensions, 
« 80% of the corresponding stellar particles are enclosed. 
In Fig. [ 9 ] we show the selected region of the (L 2 , L z ) plane 
in which global dynamics will be displayed. We also include 
the corresponding stellar particles with their actual values 
of (L 2 ,L z ). We stress that, in the following experiments, the 
orbits of the test particles are calculated using the equations 
of motion and the first variational equations for the full tri¬ 
axial NFW model. 

A global dynamical portrait of the system after an inte¬ 
gration time of 10 Gyr is shown in the left panel of Fig. |10[ 
as a contour plot of the OFLI for a grid of 450 x 450 (yield¬ 
ing a total of 77967 initial conditions). For this timescale, 
almost the whole region of angular momentum space ap¬ 
pears regular. Only a small number of invariant manifolds 
and narrow resonances are observed. The invariant man¬ 
ifolds (separatrices) show up as arcs that separate differ¬ 
ent orbital families (the large resonance domains) while the 
small (or high-order) resonances arise as channels, the cen¬ 
tres of which correspond to a chain of stable 2D resonant 
elliptic tori and the margins of which are formed by a chain 
of 2D hyperbolic (or unstable) tori. In any case, this figure 
shows that most of the phase space seems to be populated 
by regular orbits. 

Although the perturbation is not actually small enough 
for the Aq-A2 halo model, namely, si « 0.17, £2 « 0.63, 
chaos seems to be almost irrelevant after 10 Gyr of evolu¬ 


Table 4. Ensembles of 90000 initial conditions sampled uniformly 
in a box of size ~ 10 6 , the centre of which (given in the table) has 
been identified as corresponding to a chaotic orbit (recall that in 
the figures displaying the evolution of the unperturbed integrals, 
L 2 is given on a logarithmic scale). 


Ensemble 

l°gio(L 2 ) 

L z 

(0 

5.045 

275 

(h) 

6.165 

1000 

(hi) 

6.145 

25 

(iv) 

6.405 

25 


tion. The amount of chaos observed in this experiment, as 
measured in Section |4.1.3[ is ~ 12.44%. Consequently, no 
secular variation of the unperturbed integrals (L 2 , L z ) is ex¬ 
pected, since diffusion can only occur within a connected 
chaotic region of finite size. The main effect of the perturba¬ 
tion is to generate new families of orbits that, on timescale, 
are regular box and tube orbits and subfamilies of the latter. 

To look for diffusive phenomena or chaotic mixing in 
this model, we perform a second computation of the OFLI 
over a larger timescale of 250 Gyr. Although meaningless 
from a physical point of view, this experiment serves to 
unveil chaotic motion that still appears regular at 10 Gyr, 
mainly as the result of sticky orbits. We are also interested 
in the timescale on which such diffusion takes place. 

In the right panel of Fig. [10] we show an OFLI con¬ 
tour plot for a grid of 300 x 300 (34670 initial conditions). 
This reveals some hyperbolic structures, shown in red. The 
true Arnold web is mostly unveiled - it covers a consider¬ 
able domain in phase space. Indeed, the fraction of chaos 
observed in this experiment amounts to the ~ 43.37% of the 
orbits considered. Since a connected chaotic region of notice¬ 
able size exists, some secular variation of the unperturbed 
integrals (L 2 , L z ) would be expected, giving rise to fast dif¬ 
fusion (as we discuss below). For this longer timescale, the 
invariant manifold separating box from tube orbits in the 
more regular part of phase space (labelled as separatrix) is 
more clearly outlined. In addition, the high order resonant 
structure on the right, which appeared as channels crossed 
by several narrower resonances on the plot for 10 Gyr, now 
shows up as an entangled assemblage of unstable manifolds, 
leading to strongly unstable or chaotic dynamics (red com¬ 
ponent in the right panel of Fig. 101. 

The long-term diffusion of the unperturbed integrals of 
this system are determined by the topology of all its reso¬ 
nances (the right panel of Fig. |10| ) which have been detected 
efficiently by our application of the OFLI. To gain insight 
into how this diffusion operates, and to illustrate the roam of 
the unperturbed integrals, we have traced several orbits with 
initial conditions embedded in different stochastic domains 
in the (L 2 ,L Z ) plane. The wandering of the unperturbed 
integrals has been followed over 250 Gyr for ensembles of 
90000 initial conditions sampled uniformly in boxes of size 
~ 10 -6 . The centres of these boxes (listed in Table [T]) are 
identified as the most chaotic regions by the OFLI. These 
ensembles are indicated by small green rectangles in the fol¬ 
lowing figures. The equations of motion were integrated with 
a time-step of 1 Myr. Each crossing of an orbit through a 
spherical shell of radius 0.1 kpc around the Sun is depicted 
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Figure 10. OFLI contour plots for 10 (left panel) and 250 (right panel) Gyr for the Aq-A2 halo model for (xo,yo,zo) = (8,0,0), Eo — 
—204449 km 2 s —2 . The solid black lines in the colour bars indicate the values of the threshold taken to distinguish regular from chaotic 
motion. Then, warm colours indicate chaotic motion while cool colours represent regular motion. The Arnold web is mostly unveiled 
and, as we can see from the warm colours of the right panel, it covers a considerable domain in phase space. 
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Figure 11. Long-term diffusion over 250 Gyr for ensemble (i) of 
initial conditions (depicted in green) overplotted on the Arnold 
web. Diffusion proceeds along the stochastic layer separating box 
from tube orbits. 


Figure 12. Long-term diffusion over 250 Gyr for ensemble (ii) of 
initial conditions (depicted in green) overplotted on the Arnold 
web. The unperturbed integrals remain confined to a rather small 
domain, hence diffusion turns out to be mostly inefficient. 


as a black dot in the OFLI contour plot for the 250 Gyr 
interval. 

Fig. 11 shows the time evolution of ensemble (i). Dif¬ 
fusion proceeds along the stochastic layer separating box 
from tube orbits. Recall that, even in the nearly completely 
regular scenario (for an almost vanishing perturbation, i.e. 
fj, s —> 0), the orbits will remain in an exponentially narrow 
chaotic domain around the separatrix, at least for very large 
times. 

Meanwhile, when the evolution of ensemble (ii) is con¬ 
sidered (Fig. |12| ), we recognise that the jump in pseudo¬ 
integrals between the two thin parallel curves corresponds 
to the borders of a resonance. In this case, the unperturbed 
integrals remain confined to a rather small domain, so that 
diffusion turns out to be inefficient for a rather large interval. 

We stress that the diffusion we observe has some geo¬ 
metrical resemblance to the theoretical conjecture of Arnold, 


according to which diffusion proceeds through phase space 
along the chaotic layers of the full resonance web. However, 
it is clear that Arnold’s mechanism is not the way to under¬ 
stand this diffusion, since not only it is impossible to find 
any path for the orbits, but also the perturbation is not 
small enough ( jj, s < 1). Thus the confined variation of L 2 
and L 2 we find should be interpreted in terms of another 
regime, the well known overlap of resonances. Even though 
fast diffusion can occur in this scenario, clearly this is not 
the case here. 

To illustrate how diffusion progresses, Fig. [13] shows 
snapshots corresponding to 50 (left panel) and 70 Gyr (right 
panel) for a third ensemble. We notice that diffusion ad¬ 
vances along the outermost edge of the separatrix discrimi¬ 
nating box from tube orbits, near the bottom of the figure, 
and climbs to slip over the left part of the web’s upper bor¬ 
der. When a larger time interval is considered, the points are 
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Figure 13. Diffusion over 50 Gyr (left panel) and 70 Gyr (right panel) for ensemble (iii) of initial conditions (depicted in green) 
overplotted on the Arnold web as in previous figures. Diffusion advances along the outermost edge of the separatrix discriminating box 
from tube orbits. 
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Figure 14. Long-term diffusion over 120 Gyr (left panel) and 250 Gyr (right panel) for ensemble (iii) of initial conditions (depicted in 
green) overplotted on the Arnold web as in previous figures. Diffusion spreads out to cover the full width of the resonances. 


seen to spread out to cover the full width of both resonances 
as shown in Fig. 14 for 120 (left panel) and 250 Gyr (right 
panel). 


Meanwhile, for ensemble (iv) after 30 Gyr (Fig. |15| 
left panel) the unperturbed integrals wander over the small 
chaotic sea at the bottom right corner of the web, where 
an intricate overlap of resonances is observed. For a larger 
timescale (70 Gyr, right panel) the roaming of (L 2 , L z ) is 
confined to a domain of nearly the same extent (naturally 
more populated). 


Fig.[l6]also illustrates that the orbits with initial condi¬ 
tions in ensemble (iv) sweep a bounded fraction of prime in¬ 
tegral space after 120 (left panel) and 250 Gyr (right panel). 
However, the chaotic component is far from being fully con¬ 
nected, even for this large timescale, since some chaotic do¬ 
mains still remain unexplored. 


From all these experiments we clearly see that in this 
Hamiltonian representation of the Aq-A2 DM halo, diffusion 


or chaotic mixing is completely irrelevant on any realistic 
timescale. 


5 DISCUSSION AND CONCLUSIONS 

The phase space distribution of halo stars in the Neighbour¬ 
hood of the Sun potentially holds an invaluable source of 
information about the assembly history of the Milky Way. 
Stellar streams are the most direct signals of Galactic ac¬ 
cretion and their identification in the Solar Neighbourhood 
is of fundamental importance to the study of galactic dy¬ 
namics. Our capability to detect these stellar streams might 
be threatened by local chaotic mixing processes that can 
smooth out the phase space distribution function on a very 
short timescale (Section [I]). In this work we have explored 
whether chaotic mixing can play an important role in shap¬ 
ing the phase space distribution of orbits local to the present 
position of the Sun Our results reinforces the idea that this 
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Figure 15. Diffusion over 30 Gyr (left panel) and 70 Gyr (right panel) for ensemble (iv) of initial conditions (depicted in green) 
overplotted on the Arnold web as in previous figures. The unperturbed integrals wander over the small chaotic sea. 
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Figure 16. Long-term diffusion over 120 Gyr (left panel) and 250 Gyr (right panel) for ensemble (iv) of initial conditions (depicted in 
green) overplotted on the Arnold web as in previous figures. Diffusion sweeps a bounded fraction of the prime integral space —however, 
the chaotic component is far from being fully connected. 


process is very inefficient within a physically meaningful 
timescale, even within the Solar Neighbourhood. 

The degree of substructure in the Solar Neighbour¬ 
hood’s phase space distribution depends on several fac¬ 
tors. First, dynamical timescales in the inner regions of the 
Galaxy are relatively short. In general, stellar streams within 


this region are expected to be spatially well-mixed (Gould 


20031. In addition, due (at least) to the triaxial and cuspy 


nature of the underlying gravitational potential, a fraction of 
local streams could be evolving on chaotic orbits. Chaos, in 
the Lyapunov sense, indicates exponential divergence of ini¬ 
tially nearby orbits in phase space. Stream stars with chaotic 
orbits would experience very rapid mixing, with local spatial 
densities decreasing at an exponential rate such that their 
detection at the present day would be unlikely. Furthermore, 
regions filled with chaotic orbits can foster chaotic diffusion, 
which erases the ‘dynamical memory’ imprinted in phase 
space and effectively produces a smooth distribution func¬ 


We have shown that, even though all of these processes 
are undoubtedly active in the Solar Neighbourhood, they do 
not necessarily imply a significantly smooth phase space dis¬ 
tribution function. One of the key factors in this discussion 
is the relevant timescale of the system, which serves as an 
upper bound for the actual time available to develop chaos. 
In particular, the (dark matter halo) potential in the inner 
20 kpc of the Galaxy is not expected to have evolved sig¬ 
nificantly within the last ~ 8 Gyr (z = 1, see | Wang et al.| 


2011 Buist & Helmi 2014 1 . The efficiency of chaos should 


be evaluated with regard to this timeframe. 

To characterise the true impact of chaos in shaping lo¬ 
cal stellar halo phase space structure, we used fully cos¬ 
mological IV-body simulations of the formation of Milky 
Way- like dark matter haloes, coupled with a semi-analytic 


tion (Section 4.2.11. 


model of galaxy formation (Section 2.11, to sample the 
phase space distribution of Solar Neighbourhood-like vol¬ 
umes (Section |2.3| ). We modelled the dark matter halo po¬ 
tential with a triaxial extension of the well known NFW 
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density profile (Section 2.21 and used it to integrate the 
equations of motion, coupled with the first variational equa¬ 
tions, to compute the Orthogonal Fast Lyapunov Indicator, 
OFLI (Section 2.41. This chaos indicator allowed us to ro¬ 
bustly characterise the dynamical nature of stellar particles. 
Orbits were classified into three different components: reg¬ 
ular, sticky and chaotic. An important difference between 
these three families is the rate at which the local (stream) 
density around a given particle decreases as a function of 
time. While for regular orbits the local density decreases as 
a power law, for chaotic orbits it does so at an exponential 
rate. In between we find the so called sticky orbits. Sticky 
orbits behave as regular for a given period of time, after 
which their behaviour becomes chaotic. 

In all Solar Neighbourhood-like volumes analysed we 
find that, even within these strongly triaxial potentials, only 
< 20% of the stellar particles reveal their chaotic nature 
within a Hubble time. We find that ~ 30% of orbits can be 
characterised as regular. For the remaining ~ 50%, namely 
‘sticky’ orbits, it takes in general much longer than 10 Gyr 
to reveal their chaotic behaviour. The fraction of sticky or¬ 
bits is particularly important because, for halo stars moving 
on such orbits, chaotic mixing may not have enough time to 
operate even though the orbits themselves have an intrin¬ 
sically chaotic nature. It is important to mention that, in 
all cases, we have considered a simplified representation of 
the underlying galactic potentials. To explore whether this 
approximation could affect our results, orbits of the stellar 
particles associated with the Aq-A2 local volume were in¬ 
tegrated on potentials with axis ratios extracted from the 
remaining four dark matter haloes. In all cases, we kept the 
virial mass and concentration parameters fixed to those as¬ 
sociated with halo Aq-A2. If our results would be highly 
sensitive to the shape of the potential, this change to the tri- 
axiality should significantly increase the fraction of chaotic 
orbits. In all cases we find that the fraction of orbits present¬ 
ing chaotic behaviour within 10 Gyr is < 30%. The obtained 
fractions are very similar to the one obtained with the axis 
ratio extracted self-consistently from the dark matter halo 
Aq-A2. This indicates that, as long as the main sources of 
chaos are included in the model (i.e., central cusp, triaxial 
shape and its radial dependence), slight variations of the 
galactic potential should not dramatically alter the global 
dynamics of the system. 

Our results indicates that chaotic mixing, although 
non-negligible, is not a significant factor in erasing local 
signatures of accretion events. This is in agreement with the 
results of G13, who quantified the number of stellar streams 
in the same Solar Neighbourhood like volumes considered 
in this work. Their results suggest that the strongest limita¬ 
tion on quantifying substructure is mass resolution, rather 
than diffusion due to chaotic mixing. They found that, in 
the best-resolved local volumes, the number of identifiable 
streams ranges from « 300 to 600, in very good agreement 
with previous analytic predictions ( |Helmi fe White||1999[ 
Helmi et al.||2003l. It is important to note that the orbital 


classification presented in G13 most likely overestimates the 
fraction of stellar particles moving on chaotic orbits. As a 
consequence of the high but still finite particle resolution of 
these simulations, G13 were forced to track the time evolu¬ 
tion of local (stream) densities within relatively large spheres 
(R > 4 kpc). As shown in this work, such large spheres are 


likely to encompass stellar particles that exhibit very differ¬ 
ent dynamical behaviours, and hence do not faithfully repre¬ 
sent the time evolution of their corresponding local densities. 

Despite the optimistic view described above, some rel¬ 
evant caveats of the present work must be discussed and 
addressed in follow-up work. First of all, the galactic po¬ 
tentials considered in this work are associated only to the 
underlying distribution of dark matter. In a more realistic 
model, this potential should also account for the mass distri¬ 
butions of the galactic disc, bulge and a supermassive central 
black hole (e.g. Siopis fe Kandrup||2000 Kandrup & Siopis 


2003 Valluri et al.||20l0 20121. Note, however, that previ¬ 

ous studies including multi-component galactic potentials 
have successfully identified large amounts of substructure in 


Solar Neighbourhood-like volumes (e.g. Helmi & de Zeeuw 


2000 Gomez et al.|2010 Sharma fe Bland-Hawthorn|2011 1 . 
Additionally, studies based on stellar haloes obtained from 
fully cosmological hydrodynamical simulations find an over¬ 
all fraction of chaotic orbits that, in all cases, is < 20%, 


in agreement with our results (Valluri et al. 2010 20131. 


We have also considered static potentials that are not al¬ 
lowed to evolve as a function of time. While an evolving 
potential could enhance the efficiency of diffusion in phase 
space (see for instance Penarrubia|201 31, previous attempts 
to characterise the degree of substructure in local volumes, 
taking into account the variation of the Galactic potential 
in a cosmological context, have suggested that this effect 
may not be significant after all. As previously discussed, at 
least within its inner regions, the Galactic potential is not 
expected to have evolved significantly during the last ~ 8 
Gyr. We will further explore the validity of the assumptions 
adopted in this work in a follow-up study (Cincotta et al. 
2015, in preparation). 
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APPENDIX A: EXPONENTIAL DIVERGENCE 
AND VARIATIONAL CHAOS INDICATORS 

The chaos indicators, which are based on the time evolution 
of the deviation or tangent vector to the flux describing a 
given dynamical system, measure the rate of divergence of 
two initially close solutions. In order to quantify the rate 
of divergence of such two nearby orbits, let us consider a 
Hamiltonian system defined on a differentiable manifold, the 
energy surface in the present case. 

If we denote by R(p,q) the Hamiltonian with p, q 6 
K^, the energy surface is thus defined by Mh = {x : 
R(p,q) = h }. Further, on introducing the notation 


x = (P- q) ^ Mh, f(x) = (-dR/<9q, dH/dp) G Mh, 
the equations of motion can be recast as: 

x = f(x), (Al) 

so that the first variational equations take the form 


<9f 

w = —w, 
ox 


(A2) 


where w is the deviation vector and <9f/9x denotes the Ja¬ 
cobian matrix of f. Let 7(t) (an orbit) denote a solution of 
(All for the initial condition xo G Mh- Introducing some 
norm in Mh, II • II, we denote 


o(t) 


ll w MI 

ll w o|| 


which characterises the Hamiltonian flow in a small neigh¬ 
bourhood of 7 (t). Therefore the mean local rate at which 
nearby orbits to 7(f) diverge is given by the largest Lya¬ 
punov Characteristic Exponent (1LCE), defined as: 


1 1 Z 100 At 

1LCE 7 = lim - In 5 1 (t) = lirn - / —At. (A3) 

t-S-oo t t->oo t J 0 S'l 

The 1LCE allows to determine whether an orbit is regular or 
chaotic (that is, stable or unstable in the Lyapunov sense). 
Indeed, only when the S(t) increases exponentially fast with 
time the 1LCE would be different from zero. The inverse of 
the 1LCE, the Lyapunov time, provides the timescale for the 
manifestation of the local instability (that is to say the time 
needed for two nearby orbits to diverge by a factor of one 
e-folding). Whether 1LCE 7 is null or positive, 7 is said to 
be regular or chaotic, respectively. 

The numerical value of the 1LCE for a large but finite 
time T, is the so-called Lyapunov Indicator (LI), which is 
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the most widely used technique of chaos detection. Clearly 
the LI is a finite-time approximation of 1LCE 7 given by 


Eq. (A31. Therefore, a given orbit will be classified as either 
regular or chaotic whether LI converges to zero or to a pos¬ 
itive value, respectively. The inverse of the LI is the finite¬ 
time approximation to the above defined Lyapunov time (a 
detailed discussion on the theory and numerical computa¬ 
tion of the Lyapunov Characteristic Exponents and, partic¬ 
ularly of the 1LCE, can be found in the extensive review of 
Skokos|20 10). 


APPENDIX B: STICKY ORBITS 

Since sticky orbits are important in the present study, we 
briefly recap their typical behaviour. The phenomenon of 
stickiness is clearly seen in near-integrable Hamiltonian sys¬ 
tems with low-to-moderate perturbations. In this direction, 
the KAM theory (for a non-rigorous approach see for in¬ 
stance, |Chirikovj|197 9; Lichten berg fc Lieberman||1983[ ) en¬ 
sures that most of the original tori associated with the inte- 
grable system survive in the presence of a sufficiently small 
perturbation. These are the irrational tori, those that satisfy 
the so-called Diophantine condition. Meanwhile, tori close 
to a resonance condition are destroyed, leading to unstable 
chaotic motion. In systems with two degrees of freedom, the 
two-dimensional invariant tori divide the energy surface (of 
dimension 3). For systems of higher dimensionality the sce¬ 
nario is much more complicated, since, for instance the KAM 
tori no longer divide the energy surface and the stickiness 
phenomena, although still present, requires further explana¬ 
tion. Therefore, chaotic orbits cannot enter a (resonant) sta¬ 
bility domain due to the presence of an invariant curve that 
act as a barrier for chaotic motion and so remain confined to 
a finite width stochastic layer around the island. By slightly 
increasing the perturbation strength, invariant curves can be 
broken down and the barriers become only quasi-barriers to 
chaotic motion. This is due to the intricate structure of the 
former invariant curves (Aubry 1983; Mac Kay et al. |[T984 ), 
an infinite set of unconnected infinitesimal parts of the ear¬ 
lier curve, strictly speaking a cantor set or cantori (see also 


Tsiganis et al. 2000 for a qualitative description). There¬ 


fore, chaotic orbits starting in a large chaotic domain could 
in general avoid the cantori. Meanwhile, a chaotic orbit with 
initial conditions close enough to the cantori might cross 
the quasi-barriers of the cantorus and stick to the stabil¬ 
ity island associated with the resonant orbit. Hence, that 
chaotic orbit would mimic a regular one in the island. Such 
a chaotic orbit is called ‘sticky-chaotic’ or simply ‘sticky’. In 
general, sticky orbits are trapped within thin chaotic layers 
and might visit the neighbourhood of different stability is¬ 
lands for a rather long time before they escape (through the 
cantori) into the chaotic sea (in the context of Poincare sur¬ 
faces of section). A sticky orbit therefore looks like a regular 
orbit for a rather long time, until its chaotic nature is clearly 
exposed - in other words its chaos onset time is relatively 
large. 
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